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TECHNIQUES  FOR  THE  EXTRACTION  OF  WATER  DEPTH  INFORMATION 
FROM  LANDSAT  DIGITAL  DATA 

1 

INTRODUCTION 

This  report  was  prepared  under  contract  DMA80Q-77-C-0053  as  part 
of  a  continuing  program  supported  by  the  Defense  Mapping  Agency, 
Hydrographic/Topographic  Center^for  the  exploitation  of  Landsat  data 
for  updating  ocean  charting  of  the  world.  Particular  emphasis  has 
been  placed  on  developing  computer  techniques  to  best  calculate  and 
extract,  reliable  water  depth  measurements  from  NASA  supplied  digital 
data  taken  from  Landsat  sensors  over  the  shallow  seas  that  are  hazardous 
to  shipping.;  Previous  demonstrations  of  the  feasibility  of  using  Landsat 
digital  data  has  been  reported  in  the  DMA  sponsored  report  entitled 
"Demonstration  of  Satellite  Bathymetric  Mapping",  ERIM  Report  No. 

1 22200- 1 -F ,  1977.  In  order  to  extract  the  measurement  of  water  depth 
most  accurately  special  attention  must  be  paid  to  the  radiometric  and 
geometric  properties  of  the  Landsat  sensor  data  as  well  as  the  correction 
of  certain  signal  variations  in  the  6  detector  arrays  used  to  scan  the 
oceans. 

Furthermore,  in  order  to  reduce  the  effects  of  variable  radiance 
contributions  to  the  received  signal  from  atmospheric  and  oceanic  sources 
an  investigation  of  mul tispectral  and  mul ti temporal  techniques  was  begun 
during  this  period.  Finally,  a  test  of  water  depth  accuracies  as 
extracted  from  Landsat  data  using  algorithms  developed  to  date  was  made 
based  on  measurements  taken  aboard  the  R/V  Constance  in  October  1977 
in  the  Bahamian  photo-bathymetric  test  range  defined  by  John  Spinning 
and  James  Hammack  of  the  Advanced  Technology  Division,  Systems  and 
Techniques  Directorate,  Hydrographic  and  Topographic  Center,  DMA. 

The  progress  to  date  (September  1978)  on  each  of  six  tasks  is 
reported  in  tne  following  sections. 
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TASK  1  -  RADIOMETRIC  CORRECTIONS  ' 

The  six-detector  scanning  configuration  used  in  the  Landsat 
mul tispectral  scanner  subsystem  has  resulted  in  continuing  problems  due 
to  calibration  differences  among  the  detectors.  These  problems  are 
caused  mainly  by  variations  in  detector  response  characteristics,  but 
are  compounded  by  the  system  design  and  processing  procedures  used  at 
NASA/Goddard.  In  this  section,  the  discussion  will  focus  on  some  of  the 
algorithms  which  have  been  developed  for  reducing  the  striping  problems 
in  Landsat  data  sets  which  have  been  routinely  processed  and  distributed 
by  NASA.  In  the  following  section,  the  NASA  processing  procedures  will 
be  discussed,  and  recommendations  will  be  made  for  modifying  these 
procedures  or  for  implementing  alternative  procedures  which  are  more 
suitable  for  bathymetric  applications. 

Several  different  types  of  destriping  algorithms  have  been 
developed,  each  with  its  own  advantages  and  disadvantages.  The  simplest 
type  of  algorithm  is  an  offset  correction,  in  which  a  constant  c^  is 
added  to  the  signals  for  detector  i  (i  =  1,  6)  in  order  to  minimize 
the  variance  of  the  entire  data  set.  The  second  type  of  algorithm,  in 
order  of  complexity,  is  a  correction  of  the  form 


where  V'^  is  the  corrected  signal  for  detector  i,  is  the  uncorrected 
signal,  and  a^  and  b-  are  parameters  which  may  be  interpreted  as 
correction  factors  for  differences  in  gain  and  offset,  respectively, 
among  the  detectors.  This  algorithm  can,  in  principle,  correct  for 
detector  differences  over  a  wider  range  of  signals  than  the  simple 
offset  correction,  assuming  that  the  detector  responses  are  linear  with 
radiance.  A  modification  of  this  algorithm  applies  a  piecewise  linear 
correction  to  each  detector,  on  the  assumption  that  the  detector  responses 


2 


Term 


are  nonlinear  but  may  be  approximated  by  a  piecewise  linear  function. 

The  third  type  of  algorithm  is  one  which  replaces  each  signal  value  by 
a  corrected  value  using  a  lookup  table.  This  procedure  is  capable  of 
making  any  kind  of  linear  or  nonlinear  correction  and  is,  therefore, 
the  most  general  type  of  algorithm. 

An  important  criterion  in  evaluating  the  various  types  of  destriping 
algorithms  is  the  accuracy  with  which  the  correction  parameters  can  be 
determined  for  a  given  scene.  For  all  but  the  simple  offset  correction, 
the  determination  of  these  parameters  requires  that  the  scene  must 
contain  areas  of  a  uniform  or  a  slowly  varying  signal  over  a  wide  range 
of  signals.  These  conditions  do  not  occur  in  some  oceanic  scenes,  for 
example  tnose  containing  large  expanses  of  open  water  and  only  a  few 
reefs  or  shoals.  For  such  scenes,  the  offset  correction  is  the  only 
useable  method  since  the  parameters  required  for  the  other  algorithms 
cannot  be  accurately  determined.  Fortunately,  an  offset  correction  with 
parameters  determined  from  the  open  water  signal  is  usually  adequate 
for  water  depth  processing,  since  the  sensitivity  to  noise  is  greatest 
for  signal  values  near  the  deep  water  signal. 

A  second  factor  to  be  considered  in  selecting  a  destriping  algorithm 
is  the  possibility  of  variations  in  the  striping  pattern  from  one  part 
of  the  scene  to  another.  This  effect  is  illustrated  in  Figure  1,  which 
shows  mean  values  for  each  detector  over  open  water  in  two  different 
parts  of  Landsat  frame  11249-14435.  The  areas  each  contain  a  total  of 
6000  pixels  (1000  pixels  per  detector),  and  are  separated  by  720  lines. 
The  difference  in  the  striping  patterns  for  these  two  areas  is  most 
evident  in  MSS5;  in  area  A,  detectors  1  and  2  are  about  1  count  above 
the  mean,  while  in  area  B,  detector  1  is  high  and  detectors  3  and  4  are 
about  1  count  below  the  mean.  An  optimum  striping  correction  for  area  B 
actually  increases  the  variance  among  detectors  for  area  A.  In  order 
to  remove  the  striping  from  the  entire  scene,  the  scene  must  be  broken 
into  subregions  within  which  the  striping  pattern  is  uniform  within 
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FIGURE  1.  DETECTOR  MEAN  VALUES  OVER  OPEN  WATER  FOR  TWO  AREAS  IN  LANDSAT 
FRAME  11249-14435,  STRIP  3.  Area  A:  Lines  325-354,  Points  1-200.  Area  B:  Lines  1045-1074 

Points  1-200. 
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acceptable  limits,  or  a  line-by-line  correction  must  be  made.  The 
former  procedure  can  be  used  with  any  destriping  algorithm,  subject 
to  the  constraints  discussed  in  the  previous  paragraph,  while  a 
line-by-line  correction  can  be  reliably  done  only  with  the  offset 
algorithm  because  the  parameters  required  for  the  other  algorithms 
cannot  be  reliably  obtained  from  a  single  line  of  data. 

A  more  detailed  discussion  of  each  type  of  correction  algorithm 
is  presented  in  the  following  subsections. 

2.1  OFFSET  CORRECTION 

The  first  step  in  obtaining  the  offset  correction  parameters  is  to 
compute  the  mean  signal  values  for  each  detector  (i.e.,  each  sixth 
line)  over  a  selected  area.  This  area  may  include  a  range  of  signals, 
if  each  detector  views  the  same  distribution  of  radiances,  but  for 
hydrographic  applications  the  best  results  are  obtained  if  an  area  of 
clear  deep  water  is  selected.  The  rationale  for  this  selection  is 
three-fold:  first,  since  such  an  area  presents  a  uniform  radiance,  the 
requirement  of  equal  radiance  distributions  for  all  detectors  is  easily 
met  even  for  relatively  small  areas.  Second,  since  the  simple  offset 
correction  guarantees  removal  of  detector  differences  only  for  a  limited 
range  of  signals,  the  selection  of  a  deep  water  area  insures  that  an 
optimal  correction  will  be  obtained  for  signals  near  the  deep-water 
signal,  where  the  depth  extraction  algorithms  are  most  sensitive  to 
noise.  In  terms  of  the  depth  error  incurred,  a  detector-to-detector 
variation  of  several  counts  near  the  upper  limit  of  the  signal  range 
may  have  less  effect  than  a  variation  of  one  count  near  the  deep-water 
signal.  The  third  reason  for  the  selection  of  a  deep  water  area  for 
calculating  striping  coefficients  is  that  the  deep  water  signals  are 
needed  for  subsequent  depth  processing,  so  the  collection  of  statistics 
for  deep  water  serves  a  double  purpose. 

After  the  detector  mean  values  have  been  obtained,  a  set  of  offset 
corrections  (c^)  is  determined.  The  criterion  for  an  optimum  correction 
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is  that  tne  variance  among  detectors  is  minimized.  This  variance  is 
defined  oy 

=  H  «i  *  c,  ■  v)2 

°i=l  1 

where  ^ 

7  =  e  <vi 

and  V.  is  the  mean  value  for  detector  i  and  C-  is  the  offset  correction 
factor  for  detector  i.  If  real  arithmetic  is  allowed,  the  optimum 
correction  factors  are  given  by 


and  the  variance  can  be  reduced  to  zero.  However,  if  the  corrected 
data  is  to  be  stored  in  the  same  integer  form  as  the  raw  data,  only 
integer-correction  factors  are  allowed.  In  this  case,  it  can  be  shown 
that  merely  rounding  off  to  the  nearest  integer  value  is  not  necessarily 
the  optimal  solution  [1],  but  that  the  optimal  solution  can  be  obtained 
by  adding  n/6  to  C-  before  rounding  off,  where  n  is  an  integer  from 
1  to  6.  Since  the  actual  value  of  n  which  results  in  the  optimal 
correction  cannot  be  predicted  a  priori ,  each  of  these  six  possible 
sets  of  correction  factors  must  be  generated  and  tested  to  find  the 
optimal  one. 

2.  2  OFFSET  AND  GAIN  CORRECTION 

Two  methods  have  been  developed  for  making  simultaneous 
corrections  for  both  offset  and  gain  differences  among  detectors.  The 
first  metnod  requires  that  detector  mean  values  be  calculated  for  two 
(or  more)  areas,  each  of  which  has  a  different  average  signal,  and 
within  which  the  signals  are  uniform  or  slowly  varying.  If  V-j(A)  is  the 
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mean  value  for  detector  i  in  area  A,  V(A)  is  the  "correct"  value  (e.g., 
the  average  overall  detectors),  and  V ^ ( B )  and  V(B)  are  the  corresponding 
values  for  area  B,  then  the  corrected  signal  for  detector  i  can  be 
written  as 


V‘ 


=  V(A) 


+ 


V(8)  -  V(A) 

V i ( B )  -  V.(A) 


tVi 


Vi(A)] 


wnere  V-  is  the  uncorrected  signal  [2].  This  correction  assumes  that 
the  detector  responses  are  col  inear.  If  there  are  nonlinearities,  the 
procedure  can  be  adapted  by  using  the  average  signals  in  three  or  more 
areas  and  approximating  the  non-linear  curve  by  a  piecewise  linear 
function. 

The  second  method  of  obtaining  the  parameters  for  a  linear  or 
piecewise  linear  correction  of  detector-to-detector  variations  is  to  do 
a  least-squares  fit  or  a  linear  regression  between  the  signals  for  one 
detector  and  the  signals  for  each  of  the  other  detectors  in  turn.  This 
procedure  requires  the  assumption  that  the  data  set  can  be  decomposed 
into  a  set  of  N  signal  pairs  (V-n,  vjn)>  n  =  1...N,  where  V-n  is  tne 
signal  for  detector  i  and  V^n  is  the  signal  for  detector  j,  such  that 
these  signals  correspond  to  the  same  radiance.  (The  most  natural  way  to 
form  these  pairs  is  to  use  signals  from  the  same  point  number  on  tne 
same  mirror  sweep.)  This  assumption  will  be  approximately  valid  if  the 
frequency  of  spatial  variations  is  small  compared  with  the  sampling  rate. 
This  condition  is  met  in  some  oceanic  scenes,  such  as  those  containing  a 
gently  sloping  bottom  with  a  uniform  bottom  reflectance,  but  not  in  all 
scenes.  Thus,  the  scene  must  be  carefully  evaluated  before  applying 
this  method. 

Assuming  that  the  conditions  required  for  the  application  of  this 
method  are  met,  the  correction  coefficients  are  given  by 
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N 

y  (v .  -  v. )(v .  -  v . ) 

L  x  in  i  ' '  jn  j  ' 
n  =  l 

N 

in  '  V* 

n=  I 


bi  =  Vj  -  ai  Vi 


N 


where  V .  =  rr  T  V . 

'  N  nS)  ln 


and 


v.  =  i  y  V. 

J  N  J" 


This  procedure  effectively  normalizes  each  detector  i  to  a  common 
detector  j,  assuming  a  linear  relationship  exists  between  the  two 
detectors.  If  the  relationship  is  non-linear,  a  piecewise  linear  fit  may 
oe  obtained  by  partitioning  the  data  set  into  two  or  more  signal  ranges 
and  doing  a  separate  least-squares  fit  for  each  range. 

2.3  HISTOGRAM  EQUALIZATION  CORRECTION 

The  final  type  of  destriping  algorithm  considered  is  the  histogram 
equalization  method  described  by  Rosenberg  [3].  In  this  method,  cumula¬ 
tive  histograms  of  data  values  are  generated  separately  for  each 
detector,  and  each  original  data  value  is  reassigned  a  new  value  (via 
a  lookup  table)  in  order  to  force  the  histograms  for  each  detector  to 
be  the  same.  Specifically,  the  procedure  is  to  find,  by  interpolation, 
the  signal  values  (V'^n)  Tor  each  detector  (i)  corresponding  to  a  cumula¬ 
tive  frequency  of  n  percent  (n  =  0,  1,...,  100).  These  va’ues  are 
then  averaged  over  all  detectors  to  give  a  set  of  "correct"  signal 
values  (V’n)  corresponding  to  the  input  values  V'^n).  However,  since 
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these  are  in  general  non-integer  values,  a  second  interpolation  is 
required  to  find  the  "correct"  output  values  corresponding  to  the 
integer  input  values  for  each  detector. 

The  basic  assumption  involved  in  this  algorithm  is  that  each  of 
the  six  detectors  views  the  same  distribution  of  radiances  in  the 
calibration  scene.  This  is  a  somewhat  less  stringent  requirement  than 
those  for  the  other  algorithms  described  above.  However,  this 
condition  must  be  met  for  the  entire  range  of  signals  to  be  corrected, 
since  the  correction  is  valid  only  over  the  range  of  signals  encountered 
the  calibration  scene.  In  the  example  cited  previously  of  a  scene 
containing  large  expanses  of  open  water  with  only  a  few  reefs  or 
shoals,  the  conditions  required  for  the  application  of  this  method  would 
probably  not  be  met  since  the  upper  range  of  signals  would  not  be 
equally  represented  in  all  detectors.  The  primary  advantage  of  this 
method,  assuming  that  the  required  conditions  are  met,  is  that 
corrections  can  be  made  for  highly  non-linear  deviations  such  as  those 
introduced  by  the  NASA  signal  decompression  routine. 
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TASK  2  -  RADIOMETRIC  AND  GEOMETRIC  PROCESSING  PARAMETERS 

This  section  contains  a  description  of  the  procedure  used  by  NASA 
for  processing  Landsat  data  and  preparing  computer  compatible  tapes 
(CCTs).  The  information  presented  here  was  gathered  from  various 
sources,  including  the  Landsat  Data  Users  Handbook  [4],  NASA  Landsat 
Newsletters  and  Bulletins,  other  NASA  reports  [5,6],  and  personal 
communications  with  NASA  personnel.  The  descriptions  are  not  intended 
to  be  complete  or  exhaustive,  but  focus  on  the  aspects  of  the  processing 
which  are  most  crucial  for  hydrographic  applications.  Subsection  3.1 
deals  with  the  Image  Processing  Facility  (IPF)  used  for  processing 
Landsat-1  and  Landsat-2  data,  and  subsection  3.2  describes  the  modifica¬ 
tions  to  the  IPF  planned  for  handling  Landsat-3  data.  In  subsections 
3.3  and  3.4,  recommendations  are  made  for  additional  or  alternative 
procedures  for  processing  Landsat  data  intended  for  hydrographic 
applications. 

3.1  LANDSAT  IMAGE  PROCESSING  FACILITY 

Data  from  Landsat  is  received  at  one  of  three  NASA  receiving 
stations  and  recorded  on  magnetic  tapes  which  are  brought  to  NASA/GSFC 
for  processing  and  distribution.  Processing  is  done  on  the  Image 
Processing  Fac-ility  (IPF)  which  converts  the  data  to  film  imagery  and 
also  copies  the  data  for  selected  scenes  onto  computer  compatible 
tapes  (CCTs).  The  subsystems  which  generate  film  imagery  and  CCTs  are 
functionally  separate  and  operate  independently  on  the  raw  data  received 
from  the  satellite.  The  Initial  Image  Generation  Subsystem  (IIGS),  which 
produces  the  film  imagery,  applies  its  own  set  of  radiometric  and 
geometric  corrections  to  the  data  based  on  internal  calibration  and 
orbital /platform  parameters.  The  subsystems  which  produce  CCTs  apply  a 
radiometric  correction  but  no  geometric  corrections  except  for  some 
re-ordering  of  the  data  to  compensate  for  misregistration  among  spectral 


ID 


2p 


bands  and  variations  in  scan  line  length.  Since  the  process  by  which 
CCTs  are  generated  is  of  primary  interest,  this  process  will  be  described 
in  more  detail  in  the  following  paragraphs. 

Two  stages  are  involved  in  the  process  of  generating  CCTs  from 
the  raw  Landsat  Multispectral  Scanner  data.  The  first  stage  is  carried 
out  by  the  Multispectral  Scanner  Preprocessor  (MPP),  which  basically 
reformats  selected  portions  of  the  raw  video  and  calibration  data  onto 
high  density  tapes  (HDTs).  The  second  stage  is  performed  by  the  Digital 
Subsystem  (DS),  which  reads  the  high  density  tapes  from  the  MPP  as  well 
as  image  annotation  tapes  containing  orbit  and  platform  data,  makes 
radiometric  corrections  to  the  data  and  writes  out  the  modified  data  on 
computer  compatible  tapes  (CCTs). 

In  addition  to  the  radiometric  corrections,  which  will  be  described 
in  the  following  paragraphs,  two  other  operations  are  carried  out  by 
the  DS  on  the  raw  video  data.  The  first  of  these  is  the  insertion  of 
registration  fill  characters  at  the  ends  of  each  scan  line  in  order  to 
compensate  for  the  band  misregistration  in  the  raw  data  caused  by 
sequential  sampling  of  the  detectors.  Because  of  the  delay  in  sampling 
(combined  with  the  displacement  inherent  in  the  4x6  fiber  optic 
matrix),  the  first  pixel  in  MSS4  corresponds  spatially  to  the  third 
pixel  in  MSS5,  the  fifth  pixel  in  MSS6,  and  the  seventh  pixel  in  MSS7. 
Therefore,  in  order  to  bring  the  bands  into  registration  on  the  CCTs, 
six  "fill  characters"  are  inserted  at  the  beginning  of  each  line  in 
MSS4,  four  are  inserted  in  MSS5,  and  two  in  MSS6.  This  pattern  is 
reversed  at  the  end  of  each  line  in  order  to  maintain  the  same  line 
length  for  each  band.  The  second  operation  carried  out  by  the  DS  is  a 
line  length  adjustment  to  compensate  for  variations  in  the  number  of 
samples  per  scan  line  from  scene  to  scene.  This  adjustment  is  made  by 
repeating  pixels  at  regular  intervals  in  each  line.  Because  of  the  data 
format  used  for  CCTs,  the  number  of  pixels  in  each  output  line  must  be 
a  multiple  of  24.  The  number  of  original  samples  per  scan  line  is 
typically  3216  +  lb  for  Landsat-1  and  3247  +  5  for  Landsat-2.  After 
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adding  the  six  registration  fill  characters,  the  maximum  number  of 
pixels  is  3228  for  Landsat-1  and  3258  for  Landsat-2.  Thus,  the  corrected 
line  length  is  3240  for  Landsat-1  and  3264  for  landsat-2,  and  there 
are  typically  18+6  repeated  pixels  in  Landsat-1  data  and  11+5 
repeated  pixels  in  Landsat-2  data. 

The  radiometric  corrections  applied  to  MSS  data  by  the  Digital 
Subsystem  are  a  decompression  to  linearize  data  which  has  been  trans¬ 
mitted  from  the  satellite  in  compressed  mode,  and  a  radiometric 
calibration  to  compensate  for  changes  in  detector  response  character¬ 
istics  with  respect  to  the  prelaunch  calibration  and  to  equalize 
changes  among  the  six  detectors  for  each  band.  Decompression  is  done  by 
a  table  look  up  routine  which  reads  in  values  from  0  to  63  and  writes 
out  values  from  0  to  124.  Except  for  input  values  of  2  and  3,  which  are 
both  assigned  on  output  value  of  2,  this  transformation  is  single-valued. 
However,  since  only  63  output  values  are  assigned,  there  are  "gaps" 
of  one  or  two  counts  in  the  output  signal  which  may  contribute  to  the 
striping  problem  by  accentuating  differences  among  the  detectors. 

The  radiometric  calibration  is  done  separately  for  each  detector 
and  each  scan  line,  using  the  calibration  data  which  is  recorded  at  the 
end  of  each  mirror  scan.  This  calibration  data  is  obtained  by  viewing 
an  internal  light  source  through  a  varible  neutral  density  filter  during 
the  retrace  interval.  Initially,  about  900  samples  are  recorded 
across  the  calibration  wedge,  but  a  subset  of  six  of  these  samples  are 
selected  by  the  MPP  and  copied  onto  the  HDTs  for  further  processing. 

(On  Landsat-1,  only  four  useable  samples  were  selected  in  high-gain  mode, 
since  the  first  two  were  saturated.)  During  preflight  calibration  tests, 
the  signal  values  from  these  calibration  samples  were  compared  with 
signals  obtained  by  scanning  across  an  external  standard  radiance 
source,  and  a  set  of  regression  equations  were  obtained  relating  the 
gain  coefficient  (b)  and  the  offset  coefficient  (a)  to  the  calibration 
signals  for  each  detector  and  each  band.  Thus  in  principle  it  is 
possible  to  calibrate  each  scan  line  individually,  using  the  calibration 
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signals  for  that  scan  line  to  calculate  a  and  b,  and  modifying  the 
data  using  the  equation 


where  is  the  raw  data  value  (after  decompression)  in  detector  i. 

In  practice,  because  of  noise  in  the  calibration  signals,  a  modification 

'  A 

of  this  procedure  is  used  in  which  the  values  a-j  and  bi  are  "filtered" 
or  averaged  over  a  number  of  scan  lines  before  being  used  for  averaging. 

The  radiometric  calibration  process  used  for  Landsat  has 
encountered  two  unforeseen  difficulties.  First,  the  process  has  not 
completely  succeeded  in  eliminating  calibration  differences  among  the 
six  detectors.  As  a  result  of  this  partial  failure,  a  modification  of 
the  algorithm  was  developed  and  has  been  implemented  as  of  July  6,  1977. 
In  the  new  algorithm,  the  data  is  first  calibrated  using  the  procedure 
described  above,  and  then  a  second  offset  and  gain  correction  is  made 
using  correction  factors  which  are  constant  over  each  full  frame  of 
data.  The  correction  factors  are  derived  from  an  analysis  of  Landsat 
data  over  a  period  of  time,  by  a  procedure  which  has  not  yet  been 
published  by  NASA.  The  second  difficulty  was  revealed  by  an  analysis 
of  variations  in  the  calibration  data.  The  calibration  signal  should 
be  independent  of  the  external  scene  viewed  by  the  satellite,  but  it  has 
in  fact  been  found  to  be  correlated  with  the  average  scene  brightness. 
The  reason  for  this  correlation  has  not  been  conclusively  determined, 
but  it  is  tiwught  to  be  due  to  a  hysteresis  or  memory  effect  in  the 
photomultiplier  detectors.  The  implication  of  this  finding  is  that  the 
calibration  of  the  data  is  suspect  in  areas  where  there  are  large 
variations  in  scene  brightness,  for  example  in  oceanic  areas  partially 
covered  with  clouds.  This  effect,  combined  with  the  line-by-line 
calibration  procedure,  may  also  contribute  to  variations  in  the  striping 
pattern  as  discussed  in  section  2  of  this  report. 
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3.2  MODIFICATIONS  OF  THE  IPF  FOR  LANDS AT- 3 

Several  major  modifications  are  being  implemented  in  the  Image 
Processing  Facility  for  handling  Landsat-3  data  (current  estimates  are 
for  the  new  system  to  be  operational  by  the  beginning  of  1979).  Perhaps 
the  most  significant  change  for  the  present  application  is  that 
geometric  corrections  as  well  as  radiometric  corrections  will  routinely 
be  made  before  generation  of  the  CCTs.  The  standard  mode  for  the  geometric 
corrections  will  be  cubic  convolution,  although  nearest-neighbor  resampling 
will  be  available  as  an  option.  There  are  plans  to  make  these  corrections 
using  ground  control  points  for  U.S.  images,  although  these  will 
probably  not  be  available  for  the  first  year  of  operation.  Gntil  the 
ground  control  point  library  is  completed,  and  for  most-foreign  scenes, 
the  geometric  corrections  will  be  made  using  satellite  orbital  and  attitude 
parameters.  Data  will  be  resampled  onto  a  square  57  x  57  meter  grid  and 
copied  onto  HDTs  at  NASA/Goddard.  These  HDTs  will  be  sent  to  EROS  Data 
Center  for  conversion  into  CCTs  and  film  products  for  distribution  to 
the  user  community. 

3.3  RECOMMENDATIONS  FOR  RADIOMETRIC  CORRECTIONS 

It  has  become  apparent  that  while  the  standard  method  of  processing 
Landsat  data  at  NASA/Goddard  may  be  acceptable  for  most  terrestrial 
applications,  a  more  specialized  processing  procedure  would  be  of 
benefit  for  hydrographic  applications  of  Landsat  data.  These 
applications  require  that  a  maximum  amount  of  both  radiometric  and 
geometric  information  be  preserved  during  processing,  even  at  the  cost  of 
increased  data  processing  complexity  on  the  part  of  the  user.  In  this 
subsection,  recommendations  are  made  for  alternative  radiometric 
processing  methods  for  Landsat  data.  Geometric  processing  is  described 
in  the  following  subsection.  Since  these  procedures  are  presumed  to  be 
carried  out  on  raw  data,  their  implementation  would  require  that  this 
data  be  obtained  directly  from  NASA  without  going  through  the  regular 
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processing  channels,  or  that  special  provision  be  made  to  disable  some 
of  the  routine  processing  steps  and  a  special  uncorrected  product  made 
available  to  hydrographic  users. 

The  primary  disadvantage  of  the  NASA  radiometric  correction 
algorithm  is  that  the  outputs  from  both  correction  steps  (decompression 
and  calibration)  are  transmitted  in  integer  mode,  and  are,  therefore, 
subject  to  roundoff  errors  which  might  be  significant  for  hydrographic 
applications.  If  line-by-line  radiometric  corrections  are  considered 
to  be  a  necessity,  it  is  probably  inevitable  that  at  least  the  final 
output  be  converted  into  integer  form  for  recording  on  magnetic  tape. 
However,  it  is  not  clear  that  the  decompression  output  needs  to  be  in 
integer  form,  since  this  output  is  immediately  used  in  the  radiometric 
calibration  process,  which  uses  floating-point  arithmetic.  It  is  also 
not  clear  that  a  line-by-line  radiometric  calibration  is  superior  to  a 
full  frame  cal ibration*.  If  full  frame  calibration  is  used,  it  would  be 
possible  to  copy  the  raw  data  directly  to  CCTs  and  supply  users  with  a 
24  x  64  element  look-up  table  in  floating  point  form  which  would  allows 
users  to  make  a  full  radiometric  correction  (decompression  and  detector 
calibration)  without  incurring  any  further  roundoff  errors. 

If  it  can  be  demonstrated  that  a  line-by-line  calibration  is 
necessary  it  is  recommended  that  a  modified  calibration  procedure  be 
used  on  tapes  intended  for  hydrographic  applications.  This  procedure 
would  eliminate  the  filtering  or  smoothing  of  calibration  data  and  use 
only  the  calibration  data  for  a  given  scan  line  to  calibrate  that  scan 
line.  To  reduce  noise,  all  available  calibration  data  (900  samples)  except 


*Full  frame  calibration  is  done,  for  example,  by  the  Canadian  processing 
facility.  A  test  is  being  planned  to  compare  the  results  of  the  two 
processing  methods  on  a  common  data  set,  preferably  one  with  large 
variations  in  average  brightness. 
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that  which  is  saturated  in  high-gain  mode  would  be  used  for  calibration 
instead  of  the  six  samples  currently  used.  This  would  probably  be 
prohibitively  expensive  for  routine  processing,  but  could  be  done  on  a 
limited  number  of  scenes.  If  this  were  done  by  an  external  agency 
such  as  DMA,  it  would  require  access  to  the  raw  video  tapes,  since  even 
the  intermediate  HDTs  contain  only  a  subset  of  six  calibration  samples. 

It  is  also  recommended  that  after  calibration,  the  data  be  recorded  using 
all  8  bits  per  byte  instead  of  7  as  is  currently  done.  The  only  reason 
for  reserving  one  bit  as  a  flag  is  apparently  to  indicate  registration 
fill  characters,  but  since  these  always  occur  in  the  same  locations 
it  does  not  seem  necessary  to  flag  them. 

If  the  advantages  of  line-by-line  calibration  cannot  be  demonstrated, 
it  is  recommended  that  full-frame  calibration  be  implemented  using  the 
look-up  table  method  described  above.  If  NASA  is  unwilling  to  make  such 
a  change,  it  is  recommended  that  DMA  seek  access  to  raw  (uncorrected) 
data  tapes  and  perform  their  own  full-frame  calibration  by  averaging 
the  calibration  data  over  the  entire  scene,  and  perhaps  also  by 
obtaining  statistics  from  the  video  data  as  outlined  in  section  2  of 
this  report. 

3.4  RECOMMENDATIONS  FOR  GEOMETRIC  CORRECTIONS 

Since  geometric  corrections  are  not  currently  done  before  CCT 
generation,  these  corrections  may  be  applied  by  the  user  to  the  CCTs 
as  received  (although  the  first  step  in  this  process  should  be  the  removal 
of  the  repeated  pixels  inserted  by  NASA  to  equalize  line  lengths). 

Several  software  packages  for  making  geometric  corrections  to  Landsat 
data  are  already  in  the  public  domain  including  the  Digital  Image 
Rectification  System  (DIRS)  developed  by  the  Information  Extraction 
Division  at  NASA/Goddard  [6].  Resampling  for  the  purpose  of  making 
geometric  corrections  can  be  done  by  several  methods,  including  nearest- 
neighbor  and  cubic  convolution  methods.  Studies  [7,8]  have  shown  that 
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while  cubic  convolution  enhances  the  appearance  of  the  image  and  is  able 
to  correct  for  subpixel  discontinuities,  it  has  considerable  disadvantages 
for  the  present  application.  These  disadvantages  are  that  high  spatial 
frequencies  are  attenuated,  small  features  are  distorted  by  "overshoot" 
and  spreading,  and  data  values  are  altered  in  a  manner  which  is  highly 
dependent  on  the  sampling  lattice.  The  nearest-neighbor  method,  although 
it  cannot  make  corrections  to  less  than  one-half  pixel  accuracy,  preserves 
all  of  the  original  data  values  and  has  essentially  no  effect  on  spatial 
frequencies.  The  nearest-neighbor  method  is  also  the  ~nly  method  which 
can  be  applied  to  data  which  has  been  classified  or  processed  in  such  a 
manner  that  output  data  values  are  not  linearly  related  to  radiance. 

Our  recommendation  is,  therefore,  that  data  processing  should  be  carried 
out  as  far  as  possible  on  uncorrected  data,  and  that  resampling  should 
be  done  by  the  nearest-neighbor  method  only  when  necessary  for  purposes 
of  display  or  registration  with  other  data  sources.  It  should  also  be 
noted  that  much  of  the  resampling  for  display  purposes  can  be  eliminated 
by  special  design  of  the  display  equipment  (e.g.,  non-square  picture 
elements ) . 

The  new  NASA  IPF  being  built  for  Landsat-3  data  poses  some  problems 
for  hydrographic  users,  since  the  standard  mode  of  operation  will  be 
to  resample  the  data  by  a  cubic  convolution  method  onto  a  square  57  x  57 
meter  grid.  The  option  to  use  nearest-neighbor  resampling  will  exist,  and 
would  probably  be  the  more  desirable  of  the  two  alternatives,  but 
either  method  will  almost  certainly  degrade  the  useability  of  the  data 
for  hydrographic  applications  and  make  corrections  for  detector-to- 
detector  variations  more  difficult  to  accomplish.  Our  recommendation 
is,  therefore,  again  for  DMA  to  seek  access  to  raw  (uncorrected)  data 
tapes  and  carry  out  their  own  geometric  corrections  when  necessary. 

This  is  especially  true  for  foreign  scenes  for  which  NASA  will  not  have 
ground  control  points.  Rather  than  do  a  second  geometric  correction 
using  ground  control  points  specifically  acquired  by  DMA  for  this  purpose, 
it  would  be  more  advisable  to  carry  out  these  corrections  only  once,  and 
when  necessary,  on  the  raw  data. 
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4 

TASK  3  -  MULTITEMPORAL  PROCESSING  METHODS 


One  of  the  principal  advantages  of  Landsat  as  a  source  of  data  for 
hydrographic  mapping  is  its  repetitive  coverage  of  a  given  scene  at  18 
day  intervals.  This  repetition  not  only  permits  the  identification  and 
separation  of  permanent  features  from  transient  effects,  such  as  water 
quality  and  atmospheric  variations,  but  also  may  be  used  to  reduce 
dependence  on  surface  measurements  by  allowing  the  extraction  of  water 
attenuation  parameters  from  an  analysis  of  scenes  at  different  tidal 
stages.  The  methods  and  procedures  which  have-  been  developed  for  such 
multitemporal  analyses  are  described  in  this  section.  The  application  of 
these  methods  to  actual  data  is  described  in  section  6  of  this  report. 

4.1  MULTITEMPORAL  SCENE  REGISTRATION 

The  first  stage  in  a  multitemporal  processing  sequence  is  the 
preparation  of  a  composite  data  set  in  which  correspondi ng  pixels  from 
each  separate  data  set  are  brought  into  registration  with  each  other. 

This  registration  can  be  done  directly  or  through  dual  rectification. 
Direct  registration  involves  a  superposition  of  one  data  set  onto 
another  in  a  single  step,  while  dual  recti f ication  involves  a  resampling 
of  both  data  sets  onto  a  common  grid.  Since  each  resampling  process 
incurs  some  error  or  loss  of  information,  the  direct  registration  method 
is  generally  to  be  preferred.  This  method  is  also  more  computationally 
efficient  and  easier  to  apply,  since  only  relative  control  information 
is  needed,  and  this  information  can  be  obtained  in  a  semi-automatic 
process  from  the  data  itself. 

In  order  to  superimpose  one  data  set  on  another,  a  set  of  transfor¬ 
mation  equations  must  first  be  defined  which  relate  the  coordinates  of 
one  data  set  to  the  other.  Several  forms  may  be  used  for  this  transfor¬ 
mation,  including  linear,  piecewise  linear  (affine)  and  polynomial 
equations.  The  linear  and  affine  transformations  are  of  the  form 
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where  L-j ,  P-|  are  the  coordinates  (line  and  point  numbers)  for  the  first 
data  set  and  are  the  coordinates  for  the  second  data  set.  In  the 

linear  transformation,  the  constants  in  these  equations  are  obtained 
by  doing  a  least  squares  fit  over  the  entire  set  of  ground  control 
points  for  the  scene,  and  the  resulting  equations  are  assumed  to  be  valid 
for  the  entire  scene.  In  the  affine  transformation,  the  scene  is  subdivided 
into  a  set  of  triangles  by  drawing  lines  between  the  ground  control 
points,  and  a  separate  set  of  equations  is  used  for  each  triangular  region. 
The  constants  in  each  set  of  equations  are  obtained  by  requiring  an  exact 
fit  at  each  vertex  (i.e.,  a  linear  interpolation  is  done  between  the 
nearest  three  ground  control  points).  In  the  polynominal  transformation 
a  least-squares  fit  is  done  for  the  entire  scene,  as  in  the  linear  case, 
but  terms  of  order  higher  than  one  are  included  in  the  equations. 

In  the  registration  procedure  which  has  been  developed  during  this 
contract,  a  linear  transformation  is  used  for  the  following  reasons. 

In  most  coastal  and  oceanic  scenes,  the  selection  of  ground  control 
points  is  limited  to  a  relatively  small  number  of  coastal  features, 
islands,  etc.  Furthermore,  there  is  often  the  possibility  of  a  shift  in 
these  features  due  to  tidal  state  or  erosion/deposition.  Therefore, 
it  is  felt  that  a  more  accurate  overall  transformation  can  be  obtained 
by  a  least-squares  fit  which  does  not  rely  too  heavily  on  any  single 
control  point,  but  uses  a  "concensus"  among  all  available  points.  This 
rules  out  the  affine  transformation ,  which  could  introduce  errors  in  the 
vicinity  of  a  bad  control  point.  The  1  inear  transformation  is  preferred 
over  the  polynomial  because  the  number  of  control  points  is  generally 
not  large  enough  to  accurately  determine  the  higher-order  terms,  and 
the  inclusion  of  these  terms  can  cause  large  errors  when  extrapolating 
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beyond  the  region  bounded  by  the  control  points.  In  the  data  sets  which 
have  been  processed  so  far,  no  significant  nonl ineari ties  have  been 
observed  which  would  warrant  the  inclusion  of  higher  order  terms  in 
the  transformation  equations. 

Control  points  are  obtained  by  locating  common  features  in  both 
data  sets  and  recording  the  coordinates  of  these  features.  This  may 
be  done  manually,  by  inspection  of  graymaps  or  other  displays  of  the 
data,  or  it  may  be  done  using  a  semi-automatic  method.  In  this 
method,  the  two-dimensional  cross-correlation  function  is  computed  for 
the  two  data  sets  in  the  vicinity  of  each  control  point.  The  correlation 
function  is  assumed  to  be  maximized  when  the  control  points  are  properly 
aligned.  Thus,  a  point  is  chosen  in  one  data  set  and  an  automatic 
search  can  be  made  to  locate  the  corresponding  point  in  the  other 
data  set.  In  practice  this  search  is  restricted,  in  order  to  save 
computer  time,  by  indicating  the  approximate  location  of  the  control 
point  in  the  second  scene.  The  amount  of  user  interaction  in  the 
form  of  selection  of  points  and  evaluation  of  results  is  desirable. 

Once  control  points  have  been  obtained  and  regressions  have  been 
run  to  determine  the  coordinate  transformation  equations,  the  data 
from  the  second  scene  are  resampled  and  merged  with  the  first  data  set. 
This  is  currently  done  by  the  nearest-neighbor  method,  rather  than  by 
interpolation,  for  the  reasons  outlined  in  section  3.4  above.  Thus, 
for  each  pixel  in  the  master  data  set,  the  corresponding  coordinates 
are  calculated  for  the  second  data  set,  the  pixel  nearest  to  this 
location  is  selected,  and  the  two  corresponding  pixels  are  combined 
and  written  out  on  the  output  file.  For  two  Landsat  scenes  with  four 
bands  or  channels  each,  the  output  file  would  contain  eight  channels 
of  data:  the  first  four  being  the  four  bands  of  the  master  data  set 
and  the  last  four  from  the  second  data  set. 

A  modification  of  this  procedure  may  be  used  for  registration  of 
data  from  different  sources,  such  as  aircraft  scanner  data  or  digitized 
film  imagery.  This  modification  is  to  first  resample  the  data  set 
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having  the  coarser  spatial  resolution  onto  a  grid  with  the  same 
resolution  and  orientation  as  the  finer  data  set.  Then,  the  correlation 
function  technique  described  above  is  used  to  align  the  data  sets  in 
the  vicinity  of  each  control  point.  Finally,  the  second  data  set  is 
resampled  and  merged  with  the  master  data  set.  If  the  master  data  set 
is  the  one  with  the  coarser  spatial  resolution,  and  the  resolution 
difference  is  such  that  two  or  more  of  the  smaller  pixels  are  contained 
within  the  larger  ones,  the  registration  procedure  should  average  all 
of  these  pixels  together  rather  than  pick  the  single  pixel  closest  to 
the  center  of  the  master  data  set  pixel. 

4.2  IDENTIFICATION  AND  SEPARATION  OF  TRANSIENT  EFFECTS 

After  two  or  more  Landsat  data  sets  have  been  spatially  registered, 
several  types  of  multitemporal  processing  techniques  can  be  applied. 

These  include  temporal  averaging,  change  detection,  and  removal  of 
transient  effects.  The  simplest  type  of  processing  is  temporal 
averaging,  in  which  the  data  sets  are  averaged  together  on  a  pixel-by- 
pixeV-  basis  in  order  to  reduce  random  noise  in  the  data.  Because  of 
the  exponential  dependence  of  the  signal  on  the  water  depth,  as  explained 
in  section  4.3,  this  averaging  must  be  done  in  a  special  way  to  preserve 
the  form  of  this  relationship.  A  simple  arithmetic  average  of  the 
signals  from  two  or  more  scene  dates  would  no  longer  be  exponentially 
related  to  depth  except  in  the  cases  where  the  water  attenuation 
parameters  and  the  solar  elevation  is  the  same  for  each  overpass.  In 
the  general  case,  the  exponential  relationship  can  be  preserved  by 
taking  the  product  of  each  correspond i ng  data  value  after  subtracting 
the  mean  deep  water  signal.  Alternatively,  the  exponential  relationship 
can  be  transformed  into  a  linear  one  by  taking  the  logarithm  of  the 
signal  after  the  deep-water  subtraction.  Subsequent  to  this  trans¬ 
formation  the  data  may  be  scaled  to  remove  differences  in  water  attenuation 
and  sun  elevation,  and  arithmetically  averaged  in  the  usual  way. 
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Random  noise  is  characterized  by  a  lack  of  spatial  coherence  in 
the  observed  signal  differences  between  mul ti temporal  data  sets.  This 
coherence  or  lack  of  coherence  can  probably  best  be  judged  by  inspection 
of  a  map  or  display  of  the  signal  difference  after  a  logarithmic  trans¬ 
formation  and  scaling  as  described  in  the  preceeding  paragraph.  Where 
a  consistent  spatial  pattern  appears  in  such  a  display,  a  different  type 
of  processing  other  than  temporal  averaging  is  called  for.  If  the  goal 
of  the  processing  is  to  detect  long-term  changes,  for  example  topographic 
changes  due  to  erosion  or  deposition,  the  difference  map  itself  may 
represent  the  final  product.  If,  on  the  other  hand,  the  purpose  of  the 
processing  is  to  eliminate  transient  effects  such  as  water'  turbid ity 
or  atmospheric  variations,  a  third  type  of  processing  is  necessary. 

In  this  transient  removal  algorithm,  the  data  values  from  each  data  set 
are  compared,  a  decision  is  made  as  to  the  "correct"  data  value  at  a 
given  location,  and  this  data  value  is  transferred  into  the  output  file. 

Actually  several  types  of  transient  removal  algorithms  are  conceivable. 
For  the  case  of  the  bitemporal  data  set  (comprised  of  two  temporal 
acquisitions),  the  simplest  algorithm  would  merely  select  the  lower 
of  the  two  values  on  the  assumption  that  most  transient  effects  (increased 
water  turbidity,  atmospheric  haze,  etc.)  cause  an  increase  in  the  observed 
signal.  This  algorithm  selects  the  wrong  data  value  in  tne  case  of 
cloud  shadows  or  strongly  absorbing  materials  in  the  water,  however. 

Another  possibility  would  be  to  choose  the  value  nearest  to  the  preceding 
data  value,  although  this  procedure  could  also  give  incorrect  results 
in  some  cases.  Both  of  these  procedures  could  be  improved  by  applying 
a  difference  criterion  to  determine  whether  the  changes  are  due  to 
random  noise  or  transient  effects.  For  example,  the  rms  difference 
between  the  entire  data  sets  could  be  determined  during  a  screening 
pass  through  the  data.  During  processing  the  difference  between 
each  individual  pair  of  data  values  are  then  compared  with  this  rms 
difference.  If  the  difference  is  below  the  threshold,  the  data  values 
are  averaged  together  to  improve  the  signal -to-noi se  ratio,  but  if 
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the  difference  is  larger  than  the  threshold  the  "correct"  value  is 
selected  by  one  of  the  above-mentioned  criteria. 

If  three  or  more  temporal  data  sets  are  used,  a  somewhat  more 
reliable  transient  removal  can  be  effected  by  a  "voting"  procedure. 

In  this  case,  the  central  value  could  be  selected,  or  the  central  value 
and  the  value  nearest  the  central  value  could  be  averaged  together,  or 
the  set  of  values  which  differ  by  less  than  a  given  threshold  could  be 
averaged  together.  Any  of  these  algorithms  would  give  more  accurate 
results  than  is  possible  for  the  bi-temporal  case,  since  it  is  unlikely 
that  two  samples  would  be  affected  by  transients  at  the  same  location. 

4.3  EXPLOITATION  OF  TIOAL  DIFFERENCES  FOR  WATER  PARAMETER  EXTRACTION 

Differences  in  tidal  state  between  two  data  sets  covering  the  same 
geographical  area  may  be  used  to  extract  some  information  about  water 
parameters  if  these  parameters  and  the  sun  elevation  are  the  same  for 
both  dates.  For  this  purpose  the  data  sets  need  not  be  actually  regis¬ 
tered,  but  the  relationship  between  the  coordinates  of  the  two  data 
sets  must  be  known.  The  procedure  is  to  select  areas  containing  deep 
water  (where  no  bottom  return  is  present)  and  shallow  water  (where  a 
strong  bottom  return  exists).  Corresponding  areas  are  selected  in  both 
data  sets,  and  mean  signal  values  are  computed  for  each  area.  The 
difference  between  the  mean  shallow  water  signal  and  the  mean  deep  water 
signal  is  denoted  by  AV-j  and  AV^  for  data  sets  1  and  2,  respectively. 
Assuming  that  the  tidal  state  is  known  for  both  data  sets,  and  that  the 
sun  elevation  and  water  optical  properties  are  the  same  for  both  data 
sets,  some  information  about  the  water  optical  properties  may  be  inferred 
from  the  relative  magnitudes  of  AV- j  and  AV^.  This  inference  is  made  on 
the  basis  of  the  following  argument. 

The  signal  recorded  by  Landsat  over  shallow  water  may  be  written 
approximately  as: 
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V  =  vd  +  arb  e 

where  is  the  deep  water  signal,  a  is  a  parameter  depending  on  the 
atmospheric  state  and  the  solar  irradiance,  rb  is  the  bottom 
reflectance,  0‘  is  the  apparent  solar  zenith  angle  underwater,  K  is 
the  effective  water  attenuation  coefficient,  and  z  is  the  water  depth. 
For  the  first  overpass,  we  may,  therefore,  write 


where  the  brackets  indicate  an  average  value  over  the  shallow  water 
scene.  For  the  second  overpass,  each  point  in  the  scene  has  a  water 
depth  z  +  Az,  where  Az  is  the  difference  in  tidal  height  between  the 
two  overpasses.  Thus,  the  average  signal  difference  for  the  second 
data  set  may  be  written  as 

Ti7~  /  -(1  +  sec  0'9)K(z  +  Az)\ 

AV2  =  a2  <Vb  e  v  c • 


Measurements  of  AV ^  and  AV2  can  be  used  to  extract  the  value  of  K 
without  detailed  knowledge  of  the  depth  in  the  scene  i*  the  sun  angles 
are  the  same  for  both  overpasses.  If  the  sun  angles  are  not  the  same, 
K  values  cannot  be  extracted  without  knowledge  of  z  and  there  is  no 
particular  advantage  in  the  mul ti temporal  approach. 

For  the  case  of  two  data  sets  with  the  same  sun  angle  (  '^  =  '  ?) 

we  may  write 
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The  value  of  will  differ  from  1  only  if  the  atmospheric  state  is 
different  for  Ihe  two  overpasses  --  for  most  cases  this  difference  may 
be  neglected  and  the  term  involving  ^1_  way  be  dropped  from  the  equation. 
An  evaluation  of  the  atmospheric  ^1  state  can  be  made  by  comparing 

the  signal  difference  for  two  targets  (e.g.,  beach  sand  and  deep  water) 
in  both  data  sets.  If  the  signal  difference  is  not  the  same  for  both 
overpasses,  the  value 


should  be  used  in  the  equation  for  K,  where  is  tile  signal  over  a 
bright  object  in  data  set  1,  Vdl  is  the  signal  over  a  dark  object,  etc. 

4.4  MULTITEMPORAL  DEPTH  EXTRACTION  TECHNIQUES 

Assuming  that  the  data  has  been  logarithmically  transformed  as 
described  in  section  4.2,  and  either  temporally  averaged  or  edited  to 
remove  transient  effects,  the  resultant  data  values  may  be  expressed  as 

X  =  ln(ar^)  -  (1  +  sec  ■)'  )Kz 

The  constants  in  this  equation  may  be  determined  by  measurements  of 
r^  and  K  or  by  measurements  of  z  at  several  locations.  If  the  bottom 
reflectance  is  uniform  in  the  scene,  this  equation  can  be  simply 
inverted  to  yield  the  depth  at  each  point.  If  the  bottom  reflectance 
is  variable,  a  correction  can  be  made  for  these  variations  if  the  depth 
is  shallow  enough  so  that  a  bottom-ref lected  signal  is  recorded  in  MSS5. 
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The  simplest  such  correction  assumes  a  uniform  reflectance  ratio  in 
MSS4  and  MSS5.  Thus,  the  depth  is  given  by 

X5  '  X4 

_  (1  +  sec  e1)  (k5  -  k4) 

A  more  general  correction  for  bottom  color  corrections  can  be  made  by 
using  more  complicated  functions  of  and  X^  £9]. 
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TASK  4  -  ATMOSPHERIC  AND  SURFACE  REFLECTION  EFFECTS 

A  preliminary  survey  has  been  made  of  existing  techniques  for 
separating  the  contributions  of  atmospheric  scattering  and  surface 
reflectance  from  the  subsurface  reflectance  of  the  ocean.  This  task  will 
continue  into  the  next  year  with  the  examination  of  various  Landsat 
scenes  under  different  illumination  conditions. 

Basically,  there  are  two  methods  of  extracting  information  about 
atmospheric  effects  from  mul ti spectral  scanner  data  in  the  visible  and 
near-infrared  regions  of  the  spectrum.  The  first  is  the  method  of 
examining  the  discontinuities  in  radiances  at  the  edges  of  shadows 
in  the  image.  This  method  was  originally  formulated  by  Piech  and  Walker 
[10]  as  follows:  the  radiance  of  a  given  object  with  reflectance  R  is 

L  =  aR  +  B 

in  direct  sunlight,  and 
L1  =  a' R  +  S 

in  shadow,  where  a  is  proportional  to  the  atmospheric  transmittance  and 
the  total  irradiance,  a'  is  proportional  to  atmospheric  transmittance  and 
skylight  irradiance,  and  8  is  the  path  radiance  ("air  light").  These  two 
equations  may  be  combined  into  the  equation 

.  L  =  (1  +  A-)  L'  -  —i  8 

a  a 

where  6  =  a  -  a' .  Thus,  by  plotting  L  versus  L'  for  at  least  two  objects 
with  different  reflectances  the  slope  (1  +  -4-)  and  intercept  — r  .•  may  be 
determined,  from  which  a  and  —  can  be  calculated.  If  the  reflectance  of 

1  i 

one  object  is  known,  the  absolute  values  of  t  and  t'  can  also  be  determined. 
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This  technique  can  be  applied  to  ocean  scenes  if  distinct  cloud 
shadows  exist  over  at  least  two  regions  of  different  reflectance  (e.g., 
deep  water  and  shallow  water).  It  should  be  noted  that  the  term  2  in  this 
case  is  actually  the  sum  of  the  path  radiance  and  the  surface-reflectance 
radiance.  Since  the  surface  reflects  skylight  more  or  less  specularly, 
depending  on  the  sea  state,  the  surface-reflected  radiance  will  be  the 
same  in  direct  sunlight  and  in  the  shadow  (assuming  that  the  view  angles 
are  such  that  the  image  of  the  cloud  is  not  reflected  into  the  field  of 
view).  This  component  is,  therefore,  included  in  the  value  of  2  determined 
by  this  method. 

This  method  assumes  that  the  reflectance  R  is  the  same  for  skylight 
and  sunlight,  which  is  not  precisely  true  for  ocean  scenes  but  is 
probably  accurate  enough  for  the  intended  purpose.  The  technique  is 
primarily  useful  for  determining  the  ratio  (5/V)  of  direct  illumination 
to  diffuse  illumination.  Since  volume  scattering  contributes  a  negligible 
amount  to  the  radiance  observed  by  satellite  over  the  open  ocean,  the 
value  of  3  is  for  all  practical  purposes  the  same  as  the  deep-water 
signal,  which  can  usually  be  observed  directly. 

The  second  method  of  extracting  information  about  atmospheric 
effects  is  through  an  examination  of  the  signal  in  the  red  or  near- 
infrared  region  of  the  spectrum.  At  these  wavelengths  the  water 
attenuation  is  so  high  that  a  negligible  subsurface  reflectance  exists 
in  moderately  clear  and  deep  water.  Thus,  the  total  radiance  observed 
under  these  conditions  is  due  to  atmospheric  scattering.  By  comparing 
this  observed  radiance  with  the  radiance  calculated  from  an  atmospheric 
model,  the  optical  thickness  can  be  ootained.  The  model  can  then  be 
used  to  calculate  the  path  radiance  and  transmittance  in  the  visible 
region  in  order  to  make  corrections  for  atmospheric  variations  in  the 
scene.  Models  and  procedures  have  been  proposed  by  several  authors 
L 1 1 ,  12,  13],  and  the  concensus  of  opinion  appears  to  be  that  MSS6  is 
the  most  useful  band  for  extracting  atmospheric  parameters.  These 
procedures  will  be  evaluated  during  the  next  year,  and  an  attempt  will 
be  made  to  develop  a  method  of  automatically  correcting  for  atmospheric 
variations  in  oceanic  scenes. 
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TASK  5  -  MULTITEMPORAL  PROCESSING  OF  BAHAMA  DATA 

Ihe  groundwork  has  been  laid  for  mul ti temporal  bathymetric  proces¬ 
sing  of  Landsat  data,  and  plans  have  been  made  to  process  several  data 
sets  from  orbit  number  232  covering  part  of  the  Great  Bahama  Bank. 

However,  because  of  delays  in  the  delivery  of  these  data  sets,  all  of 
the  anticipated  processing  has  not  been  completed  at  the  time  of  writing 
this  report.  Processing  will  continue  as  these  tapes  are  received,  and 
all  products  and  reports  of  this  processing  will  be  delivered  upon 
completion  of  this  task. 

In  the  interim,  some  mul ti tempora 1  processing  has  been  done  for  two 
frames  of  data  for  this  area.  The  two  frames  selected  for  this  prelimi¬ 
nary  study  are  a  low  gain  data  set  collected  on  29  December  1974  (frame 
10889-15033)  and  a  high  gain  data  set  collected  on  24  December  1975 
(frame  11249-14435).  Because  of  the  similarity  in  the  dates  the  sun 
angle  are  nearly  the  same  for  both  data  sets.  The  first  data  set  was 
tarken  about  two  hours  after  high  tide,  and  the  second  was  about  two 
hours  before  high  tide,  so  the  tidal  levels  are  nearly  the  same  for  the 
two  data  sets  also. 

Ihe  first  step  in  the  processing  of  these  data  sets  was  the  extraction 
of  relative  control  points  and  registration  of  the  two  scenes  using  the 
technigues  described  in  section  4.1  of  this  report.  The  next  step  was 
a  scaling  of  the  low  gain  bands  to  make  them  commensurate  with  the  high 
gain  data.  Figures  2-5  show  plots  of  the  low  gain  and  high  gain  MSS-4 
and  MSS-5  data  along  two  lines  near  North  Cat  Cay  and  South  Cat  Cay  in 
the  Bahamas.  Figures  6  and  7  show  two-dimensional  histograms  of  data 
values  in  the  low  gain  bands  versus  those  in  the  high  gain  bands.  A 
least-squares  fit  was  done  between  the  low  and  high  gain  data,  resulting 
in  the  following  equations: 
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HIGH  GAIN  MSS  4  =  0.93  +  2.32  (LOW  GAIN  MSS  4) 

HIGH  GAIN  MSS  5  -  3.46  +  3.30  (LOW  GAIN  MSS  5) 

ihis  analysis  shows  a  considerable  departure  from  the  nominal  3x 
gain  difference  for  MSS4.  The  standard  deviation  of  the  deep  water 
signals  was  also  calculated  for  the  low  gain  and  high  gain  bands.  These 
figures,  shown  in  Table  1,  indicate  that  there  is  a  slightly  higher 
signa 1 -to-noise  ratio  for  the  high  gain  data  than  the  low  gain  data, 
by  a  factor  of  about  1.4  in  both  bands.  The  low  gain  data  channels 
were  scaled  using  the  above  equations,  and  the  resulting  data  were  again 
plotted  along  the  same  lines  as  in  Figures  2-5,  and  are  shown  in 
Figures  8-11. 

After  registration  and  scaling,  the  two  data  sets  were  combined 
together  using  a  transient  suppression  algorithm  to  remove  clouds  and 
whitings,  and  improve  the  s i gnal -to-noi se  ratio  of  the  data.  This 
algorithm  tests  the  difference  between  the  data  values  in  the  two  data 
sets  on  a  pixel-by-pixel  basis.  If  this  difference  is  larger  than  a 
given  threshold  in  MSS5,  a  cloud  or  whiting  is  assumed  to  be  present 
and  the  data  set  having  the  lower  data  values  is  assumed  to  be  the 
correct  one.  A  test  is  also  made  to  determine  if  the  high  gain  data  is 
saturated,  and  if  so  the  low  gain  data  is  used.  If  no  clouds  or  whitings 
are  present  in  either  data  set,  and  the  high  gain  data  is  not  saturated, 
the  data  values  from  both  data  sets  are  averaged  together  to  reduce  the 
random  noise  in  the  data. 

The  results  of  the  mul ti temporal  transient  suppression  algorithm 
are  shown  in  Figures  12  and  13.  Figure  12  shows  the  raw  MSS-4  data 
for  the  two  frames  and  the  processed  output  on  the  right.  Figure  13 
shows  the  raw  and  processed  data  for  MSS-5.  Some  residual  effects  appear 
to  be  present  in  areas  covered  by  clouds  on  one  of  the  data  sets.  These 
effects  are  due  to  cloud  shadows,  which  are  not  removed  by  this  aloo- 
rithm,  and  to  the  inclusion  of  marginally  clouded  pixels  which  fall 
within  the  threshold  for  cloud  detection. 


TABLE  1 


DEEP  WATER  STATISTICS  FOR  MULTITEMPORAL 
DATA  SET  (lines  1730-1750,  points  376-426) 


Frame  1889-15033 
ii'iSS4:  mean  45.78 

Std.dev.  2.45 


Frame  5249-14435 
47.99 
1.60 


Combined  Data 
47.16 
1.47 


MSS5:  mean  26.91 


26.82 


Std.  dev.  1.84 


1.32 


27.02 

1.24 


MSS6:  mean  3.77 

Std.  dev.  0.57 


4.20 

0.71 


4.21 

0.50 


MSS7 :  mean  0.40 

Std.  dev.  0.51 


0.46 

0.54 


0.59 

0.53 
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of  Scaled  Low  Gain  and  High  Gain  MSS4  Data  North  of  North  Cat  Cay, 
Bahama  Bank. 
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Lay,  Great  Bahama  Bank 


(a)  Frame  10889-15033  <b)  Frame  11249-14435 

(Low  Gain)  (High  Gam) 


if)  Combined  Frames, 
with  Transient  Suppressi 


FIGURE  1 


2.  B&W  IMAGE  DISPLAYS  OF  RAW  AND  PROCESSED  MSS-4  DATA  V- 
FRAMES  10889-15033  (STRIP  3)  AND  1  1 249 -  1 -4433  (STRIP  3>. 


^Trjm 


(a)  Frame  10889-15033 
( Low  Gain) 


(b)  Frame  11249-14435 
(High  Gain) 


(c)  Combined  Frames 
with  Transient  Suppression 


FIGURE  13.  B&W  IMAGE  DISPLAYS  OF  RAW  AND  PROCESSED  MSS-5  DATA  FOR  FRAMES 
10889-15033  (STRIP  3)  AND  11249-14435  (STRIP  3). 
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Plots  of  the  mul titemporally  processed  MSS4  and  MSS5  data  along 
the  same  transects  north  of  North  Cat  Cay  and  south  of  South  Cat  Cay  are 
shown  in  Figures  14  and  15.  A  visual  comparison  with  Figures  8-11  shows 
an  apparent  reduction  in  noise,  which  is  confirmed  by  statistics  taken 
over  deep  water.  These  statistics  are  summarized  in  Table  1.  The 
standard  deviation  of  the  deep  water  signals  for  the  processed  data 
set  is  40  percent  lower  than  that  for  the  scaled  low  gain  data  set,  and 
8  percent  lower  than  that  for  the  high  gain  data  set  in  MSS4.  In  MSSb, 
this  improvement  is  33  percent  over  the  low  gain  and  6  percent  over  the 
high  gain  data  sets. 

dsing  the  criterion  that  the  maximum  penetration  depth  is  the 
depth  at  which  the  bottom  reflected  signal  is  equal  to  one  standard 
deviation,  the  penetration  depths  for  MSS4  are  21.9  meters  for  the  low- 
gain  data  set,  24.8  meters  for  the  high  gain  data,  and  25.4  meters  for 
the  multitemporal  data  set.  The  corresponding  figures  for  MSS5  are 
5.6,  6.1,  and  6.2  meters,  respectively.  Deeper  penetration  can  be 
obtained  by  spatial  filtering,  but  at  the  cost  of  decreasing  the  spatial 
resolution.  The  increased  accuracy  and  penetration  depth  obtained  by 
mul ti temporal  processing  is,  of  course,  accompanied  by  the  benefit  of 
obtaining  a  cloud-free  image. 
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Figure  14.  Plots  of  Mul t i tempora 1 1  y  Processed  MSS4  and  MSS5  Data  North  of  North 
Cat  Cay,  Great  Bahama  Bank 


7 

FIELD  VERIFICATION  IN  THE  BAHAMIAN  PHOTOBATH YMETR I C  CALIBRATION  AREA 

In  order  to  test  the  accuracy  of  the  depth  extraction  methods  and 
to  determine  the  effect  of  bottom  reflectance  variations  on  the  depth 
calculations,  a  field  trip  was  planned  and  carried  out  during  October 
1977  in  the  northwestern  part  of  the  Great  Bahama  Bank.  The  geographic 
coordinates  (latitude  and  longitude)  were  determined  from  satellite 
navigation  equipment  aboard  ship  at  14  achorage  locations.  Water  depths 
were  measured  by  lead  lines  and  fathometer  soundings  at  each  of  these 
locations,  and  fathometer  transects  were  taken  between  several  of  the 
anchorage  locations.  In  addition,  bottom  photographs  were  taken  at  IS 
locations  with  filters  corresponding  to  Landsat  bands  MSS-4  and  MSS-5 
and  with  a  calibrated  reflectance  standard  included  in  each  frame.  Of 
these  18  bottom  scenes,  10  were  taken  at  anchorage  locations  with  precise 
location  information  and  8  were  taken  at  some  distance  from  the  ship. 

Thus,  a  complete  set  of  observations  (location,  depth,  and  bottom 
reflectance)  were  made  at  10  stations  in  the  test  range.  The  station 
numbers,  geographic  coordinates,  and  depths  for  these  locations  are  shown 
in  Table  2. 

Bottom  reflectances  were  obtained  from  the  photographs  by  measuring 
the  density  of  the  film  emulsion  over  each  of  the  three  known  reflectance 
panels  and  using  these  to  plot  a  curve  of  film  density  versus  reflectance. 
An  average  film  density  was  then  measured  for  the  bottom,  and  the  cor¬ 
responding  reflectance  was  obtained  from  the  calibration  curve.  The 
bottom  reflectances  at  each  station  in  the  MSS-4  filters  are  shown  in 
Table  3. 

In  order  to  compare  the  measured  depths  with  those  calculated  from 
Landsat  data,  a  set  of  equations  relating  the  Landsat  coordinates  to  utm 
coordinates  was  developed  using  ground  control  points  supplied  by  DMA. 

The  geographic  coordinates  for  each  station  were  then  converted  to  utm 
coordinates  using  the  USGS  Transverse  Mercator  Transformation  Program 
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TABLE  2 


LOCATION  AND  DEPTH  OF 
TEN  STATIONS  IN  BAHAMA  BANK  TEST  RANGE 


STATION  NUMBER 

LATITUDE 

LONGITUDE 

DEPTH 

C-5 

25°  32. 20' N 

79°  16.92'W 

9.8  m 

D-7 

25°  43.94'N 

79°  18 . 20' W 

9.1  m 

F-ll 

25°  56. 60' N 

78°  59.47'W 

9.8  m 

G-12 

25°  56. 36' N 

78°  57.77'W 

10.4  m 

H-l  3 

25°  40. 26' N 

78°  40. 76' W 

4.9  m 

1-14 

25°  40. 46 ‘ N 

76°  42.22'W 

6.7  m 

J  - 1 6 

25°  50. 48 1 N 

78°  43.82 ' W 

12.5  m 

A- 17 

25°  48.60'N 

79°  00. 50 ‘ W 

6.1  in 

B-18 

26°  01.64'N 

79°  05.  90'  W 

10.7  m 

D-  21 

25°  46.15'N 

79°  16. 73' W 

6.  1  m 
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TABLE  3 


MEASURED  BOTTOM 

REFLECTANCES  AND 

LANDSAT 

SIGNALS  IN  BAND  MSS-4 

FOR  TEN  STATIONS 

IN  BAHAMA  TEST  RANGE 

STATION 

BOTTOM 

MSS-4  SIGNAL 

MSS-4  SIGNAL 

NUMBER 

REFLECTANCE 

FRAME  10889-15033  FRAME  11249-14435 

C-5 

0.20 

24 

67 

D-7 

U.  24 

22 

63 

F-n 

0.23 

22 

58 

6-12 

0.19 

23 

58 

H-13 

0.24 

33 

86 

1-14 

U.10 

25 

65 

J-lb 

0.10 

21 

52 

A- 1 7 

0.23 

26 

70 

8-18 

0.15 

21 

53 

D-21 

0.1b 

23 

63 
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No.  W5377,  and  these  coordinates  were  converted  to  satellite  coordinates 
using  the  equations  described  above.  This  process  was  repeated  for  the 
two  Landsat  frames  mentioned  in  Section  6  of  this  report.  The  utm  to 
satellite  coordinate  transformation  equations  for  these  two  frames  are 
as  follows: 

Frame  10889-15033  (Strip  3): 

LINE  =  37932  -  .002051  (utme)  -  .01233  (utmn) 

POINT =  -181  +  .01694  (utme)  -  .003811  (utmn) 

Frame  11249-14435  (Strip  3): 

LINE  =  37882  -  .001995  (utme)  -  .01230  (utmn) 

POINT  =  -252  +  .01/09  (utme)  -  .003866  (utmn) 

The  MSS4  signals  at  the  10  stations  were  extracted  from  both  Landsat 
frames  using  this  procedure,  and  are  listed  in  Table  3.  MSS5  signals 
were  also  extracted,  but  are  not  listed  because  they  were  in  most  cases 
not  significantly  above  the  deep  water  signal. 

Next,  water  depths  were  calculated  from  the  Landsat  MSS4  signals 
for  each  station,  using  the  procedure  described  in  the  NASA/Cousteau 
report  [14].  The  water  attenuation  coefficient  was  assumed  to  be  that 
measured  during  the  Cousteau  experiment  at  the  Great  Isaac  Station: 
i.e.  =  0.0748  m  \  (Water  attenuation  measurements  were  also  made 
during  the  DMA  field  trip,  but  because  of  high  sea  state  at  the  time 
of  the  observations,  the  measurements  were  not  considered  to  be  repre¬ 
sentative  of  normal  conditions.)  A  first  order  calculation  was  made 
using  an  average  (wet)  bottom  reflectance  value  of  22  percent  for  all 
the  stations.  This  is  the  same  calibration  used  previously  for  processing 
Frame  11249-14435  [15].  The  results  of  these  calcul ations  are  shown  in 
fable  4.  Next,  a  second  order  calculation  was  made  using  the  bottom 
reflectance  values  measured  for  each  station.  The  results  of  this 
calculation  are  shown  in  Table  5.  The  equations  used  for  these  calibrations 
are  as  fol lows : 

o0 


2p 


0  s 

where  K  =  0.0748  m""' 

V  -V  =  104  r.  for  frame  1 0889- 1 5033 
os  b 

293  r^  for  frame  11249-14435 

=  16.5  for  frame  10889-15033 
46.5  for  frame  11249-14435 


The  second  order  calculation  may  be  viewed  as  an  additive  correction  to 
the  first  order  calculation,  of  the  form 

AZ  =  J_  ln/V 

2K  n  r 

b 

where  r^  is  the  average  bottom  reflectance  (0.22)  used  in  the  first  order 
calculation. 


I  he  first  order  depth  calculation  miscalculates  the  depth  when  the 
bottom  reflectance  is  different  from  the  average  value.  However,  tne 
second  order  calculation  using  the  measured  reflectances  over-corrects 
for  these  differences  in  most  cases.  The  reason  for  this  is  probably 
two-fold:  first,  at  those  stations  having  a  partially  vegetated  bottom 
(1-14  and  J- 16)  the  bottom  photogrpahs  may  have  tended  to  favor  the 
vegetation  rather  than  the  sand  background,  thus  yielding  a  bottom 
reflectance  lower  than  the  average  value  over  the  instantaneous  field 
of  view  of  the  Landsat  sensor.  Secondly,  the  effect  of  scattering  in 
the  water  is  to  reduce  the  contrast  between  dark  areas  and  surrounding 
lighter  areas,  so  that  the  effective  bottom  reflectance  is  intermediate 
between  the  actual  value  within  the  field  of  view  and  those  in  the 
surrounding  areas.  As  a  result,  the  r.m.s.  error  for  the  second  order 
calculations  (2.3  m)  is  actually  larger  than  that  for  the  first  order 
calculations  (1.6  m) .  The  r.m.s.  error  for  the  first  order  corrections 
is  about  18  percent  of  the  average  depth  for  the  ten  stations. 
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TABLE  4 

FIRST  ORDER  DEPTH  CALCULATIONS  FOR  TEN  STATIONS 


STATION 


NUMBER 

FRAME  10889-15033 

FRAME  11249-14435 

AVERAGE 

C-5 

7.  d  m 

7.7  m 

7.6m 

D-7 

9.5  m 

9.1  m 

9.3  m 

F-ll 

9.5  m 

11.5m 

10.5m 

G-12 

8.4  m 

11.5m 

10.0  m 

H-13 

2.2  m 

3.3  m 

2.8  m 

1-14 

6.6  m 

8.3  m 

7.4  m 

J-lb 

10.9  m 

16.5  m 

13.7  m 

A-l  7 

5.9  m 

6.7  m 

6.3  m 

B- 18 

10.9  m 

15.3  m 

13-.1  m 

D-21 

8.4  m 

9.1  in 

8.8m 

Term 


TABLE  5 


SECOND  ORDER  DEPTH  CALCULATIONS  FOR  TEN  STATIONS 


STATION 

NUMBER  FRAME  10889-15033 


C-5 

6.9  m 

D-7 

10.1  m 

F-ll 

9.8  m 

G-12 

7.4  m 

H-13 

2.8  m 

1-14 

1.3  m 

J-lb 

5.6  m 

A- 17 

6.2  m 

B-18 

8.3  m 

0-21 

6.3  m 

FRAME  11249-14435 

AVERAGE 

7.1  m 

7.0  m 

9.7  m 

9. 9  in 

11.8  m 

10.8  m 

10.5  m 

9.0  m 

3.9  m 

3.4  m 

3.0  m 

2.2  m 

11.2  m 

8.4  m 

7.0  m 

6.6  m 

12.7  m 

10.5  m 

7.0  m 

6.6  m 

d3 


7.1  EFFECTIVE  REFLECTANCE  DERIVED  FROM  LANDSAT 


To  gain  further  insight  into  the  problem  of  varying  bottom  reflect¬ 
ance,  the  effective  reflectance  was  calculated  for  each  of  tne  lu  sites 
by  assuming  the  depth  was  known  and  solving  the  depth  equation  for  r^. 

Table  6  shows  these  calculations  for  both  scenes. 

For  the  frame  10889-15033,  the  measured  reflectance  as  derived  from 
the  underwater  photograph  of  the  calibration  panel  agreed  with  the  calcu¬ 
lated  effective  reflectance  for  the  Landsat  pixel  for  5  of  the  lu 
stations.  The  remaining  cases  suggested  that  at  those  test  sites  the 
measured  reflectance  was  biased  by  too  much  vegetation  when  compared  to 
the  size  of  the  Landsat  pixel  and  consequently  lower  bottom  reflectance 
was  reported  so  that  in  the  second  order  calculation  higher  rms  errors 
were  introduced. 

A  similar  pattern  exists  for  frame  11249-14435  at  5U  percent  of 
the  stations,  the  locally  measured  reflectance  approximates  the  calculated 
effective  reflectance  within  a  15  to  20  percent  error  range.  The  use 
of  two  spectral  channel  depth  calculation  algorithms  has  been  shown  to 
be  effective  in  adjusting  for  these  variations  in  bottom  reflectance  [15]. 
Further  tests  of  the  effectiveness  of  bottom  reflectance  on  the  accuracy 
of  depth  measurements  are  planned  using  the  data  collected  in  the  summer 
of  1978  in  the  Bahamas  Photobathymetric  Calibration  range. 

A  number  of  data  sources  including  different  scale  aerial  photographs, 
mul tispectral  active/passive  scanner  data,  and  multidate  Landsat  data 
will  be  used  to  investigate  the  variable  spatial  and  spectral  aspect  of 
this  parameter. 
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TABLE  6 

COMPARISON  OF  MhASuRED  SUBPIXEL  BOTTOM  REFLECTANCE 


AND  CALCULATED 

EFFECTIVE  REFLECTANCE  FOR 

A 

LANDSAT 

PIXEL 

Station 

Frame 

Frame 

Sub-Pi xe 

Number 

I0889-Ib033(rb) 

11249-14435(rb) 

Measured 

C-8 

.22 

.30 

.20 

0-7 

.21 

.22 

.24 

F-l  1 

.23 

.17 

.23 

6-12 

.30 

.19 

.19 

H-13 

.33 

.28 

.24 

I-U 

.22 

.17 

.10 

J-16 

.28 

.12 

.10 

A-17 

.23 

.20 

.23 

B-18 

.21 

.1  1 

.15 

D-21 

.16 

.14 

.16 

5o 
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